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We develop a self-consistent finite element method to study spontaneous emission at nanoscale 
proximity of plasmonic waveguides. In the model, it is assumed that only one guided mode is 
dominatingly excited by the quantum emitter. With such one dominating mode assumption, the 

■ cross section of the plasmonic waveguide can be arbitrary. We apply our numerical method to 
' calculate the coupling of a quantum emitter to a cylindrical nanowire and a rectangular waveguide, 
, and compare the cylindrical nanowire to previous work valid in quasistatic approximation. The 

^N) ■ fraction of the energy coupled to the plasmonic mode can be calculated exactly, which can be used 

to determine the single optical plasmon generation efficiency for a quantum emitter. For a gold 
nanowire we observe agreement with the quasistatic approximation for radii below 20 nm, but for 

■ larger radii the total decay rate is up to 10 times larger. For the rectangular waveguide we estimate 
an optimized value for the spontaneous emission factor B of up to 80%. 

CN \ PACS numbers: 42.50.Pq, 73.20.Mf, 78.55.-m 
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It has long been realized that the spontaneous emission rate is not an intrinsic property of a quantum emitter itselft. 
The general explanation is that the spontaneous emission rate depends on the transition strength between the upper 
and lower level of the quantum emitter and the local density of optical states. The local density of states measures 
the available number of electromagnetic modes into which the photons can be emitted at a specific location of the 
emitter, and can be manipulated by tailoring the photonic environment of the emitter. A number of structures such 
r"| ■ as interfaces^, cavities^, photonic crystals^ and waveguides^ have already been used to modify the spontaneous 
rS i' emission rate. Apart from fundamental studies, engineering the spontaneous emission rate of a quantum emitter can 
lead to new possibilities to boost the efficiency of optoelectronic devices, i.e., single-photon sources, low threshold 
lasers, and LED-lightening. 

^ . As an alternative to dielectric materials, the spontaneous emission rate can be manipulated by subwavelength 
' metallic systems, which support surface plasmon polaritons. Surface plasmon polaritons are electromagnetic exci- 
te I tations associated with charge density waves on the surface of a conducting object. The tight confinement of the 
OA ' electromagnetic field to the metal-dielectric interface due to the boundary condition constraints gives the possibility 
\ of inventing new ways to enhance light-matter interaction, such as efficient single optical plasmon generatioitiSili, 
^\ single molecule detection with surface-enhanced Raman scatterin g^^'^'^ , enhanced photoluminescence from quantum 

_ wellsi^, and nanoantenna modified spontaneous emission— li^iil. 
0^ . Although limited by the intrinsic losses of the metals in the optical frequency range, different metallic structures 
have been extensively studied in the last few years due to the possibilities of integration and miniaturization. The 
dramatic enhancement of the field intensity due to the field concentration and geometric slowing down of the mode 
propagation provides an excellent platform to study single photon nonlinear optics^, and light matter interaction at 
^ , the single-emitter-single-photon level. There are also considerable interests in surface plasmons for in sub-wavelength 
opticai^ and applications in sensing, near field imaging, waveguiding and switching below the diffraction limit2Si2ii2^i2^. 
The study of plasmonic effects to enhance light-matter interaction and the preferential emission of, e.g., a quantum dot 
into a desired mode is currently a hot research topic. It is important for solid-state quantum information devices as 
well as for improving our understanding of light-matter interaction at the nanoscale. So far, there are a few theoretical 
papersiii^lon this topic, and they employ simplifying assumptions that limit their applicability for analysing realistic 
structures, e.g. by assuming geometrical shapes that are inconsistent with current fabrication technology and making 
assumptions that are only valid at some length scales. The realistic description of all competing, radiative and non- 
radiative, decay channels for an emitter placed in close proximity to a plasmonic waveguide is important in order to 
understand the physics and the fundamental limitations. 

The present paper focuses on modeling of the spontaneous emission of a quantum emitter at nanoscale proximity 
to the realistic plasmonic waveguides by using a finite element method, with special emphasis on calculating the 
spontaneous emission j3 factor. The p factor describes the fraction of the emitted energy that is coupled to the 
plasmonic mode. Subwavelength waveguiding of plasmons in metallic structures has been studied theoretically22,i^ 
and has also been observed in a number of recent experiments^^. Enhanced spontaneous emission of an emitter coupled 
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FIG. 1: Different emission channels involved in the decay process of a quantum emitter (red dot) coupled to a plasmonic 
waveguide. In the radiation channel the photons are traveling in free space. In the plasmonic channel the plasmonic modes are 
excited and guided by the metallic nanowire. In the non-radiative channel, electron hole pairs are generated. 

to plasmonic waveguides has been propose d^^'^^ and experimentally demonstratediS recently. Chang et al.— studied 
the spontaneous emission of an emitter coupled to a metallic nanowire by using the quasistatic approximation. Jun 
et al.^^ employed a FDTD numerical method to study the different spontaneous emission decay rates of an emitter 
coupled to a metallic slot waveguide, but with assumptions for the local density of states of the plasmonic mode. A 
self-consistent model with rigorous treatment of all the spontaneous decay rates involved, i.e. radiative as well as 
non-radiative, has not been presented in the literature. The aim of this paper is to provide such a detailed modeling. 

As shown in Fig. [U we consider an ideal quantum emitter coupled to a plasmonic waveguide. The excitation energy 
of the quantum emitter can be dissipated either radiatively or non-radiatively. Radiative relaxation is associated with 
the emission of a photon, whereas non-radiative relaxation can be various pathways such as coupling to vibrations, 
resistive heating of the environment, or quenching by other quantum emitters. The resistive heating of the metallic 
waveguide is then the only mechanism of non-radiative relaxation considered in our model. The quantum emitter is 
positioned in the vicinity of the metallic nanowire, thus there are three channels for the quantum emitter to decay 
into, i.e., the radiative channel, the plasmonic channel and the non-radiative channel. The corresponding decay rates 
are denoted by 7rad, 7p;, and ^nonrad, respectively. The radiative channel is the spontaneous emission in the form of 
far field radiation. The plasmonic channel is the excitation of the plasmonic mode, which is guided by the plasmonic 
waveguide. The non-radiative channel is associated with the resistive heating of the lossy metals, which is due to 
electron- hole generation inside the metals. The spontaneous emission (3 factor is defined by /3 = ^J''^' ^ , where "ftotai is 
the sum of the three rates, jtotai — Irad + Inonrad + Ipi- The (3 factor gives the probability of the photon "coupling" 
to the plasmonic mode, when the single quantum emitter decays. 

This paper is organized as follows. In Sec. II, the computational principle and the numerical method are presented. 
First we study the dispersion relation and the mode properties of the plasmonic waveguide, and then we calculate 
the decay rate into the plasmonic channel in a 2D model by taking advantage of the translation symmetry of the 
waveguides. Finally, the wave equation with a current source in a 3D model is solved numerically, and the total decay 
rate of the quantum emitter is extracted by calculating the normalized total power emission of the current source. 
Section III presents the results and discussion obtained by applying the numerical method to two different plasmonic 
waveguides. Section IV concludes the paper. 



II. COMPUTATIONAL APPROACH 



A. Dispersion relation and decay rate into the plasmonic channel 



The starting point of the numerical analysis of the waveguide is the wave equation for the electric field. 
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FIG. 2; Dispersion relation versus radius for the metallic nanowire . Inset (a) shows the waveguide structure. Inset (b-g) show 
field orientation of the possible eigenmodes supported by the waveguide. 




FIG. 3; Dispersion relation versus side length of the square plasmonic waveguides . Inset (a) shows the waveguide structure. 
Inset (b-e) shows field orientation of the possible eigenmodes supported by square plasmonic waveguides. 



where fco = ^^/£qJM) is the vacuum wave number, Sr denotes the relative dielectric constant and Hr represents the 
relative permeability constant, which is a constant in our model. Due to the invariance along the Z axis, the Z- 
dependence of the solution to the wave equation must be that of a plane wave (complex exponential). 

Through out the paper, j denotes For the guided plasmonic modes, at a specific frequency cu two quantization 

indices are needed to specify a complete set of orthogonal modes, i.e., a = {p, /?}. f3 denotes the propagation constant 
(the component of the wave vector along the Z-axis), and the index p represents the polarization of the mode. The 
waveguide structure examined consists of two regions f2 and A. f2 is the lossy metal core, which is surrounded by an 
infinite lossless dielectric medium A. The transverse component of the wave vector fulfills (jfci±)^ + = with 
i G [fl, A], where ki± and Si are the transverse component of the wave vector and the relative permittivity. 

The finite element method can be utilized as a numerical tool to calculate the guided plasmonic modes. The 
infinite dielectric medium is truncated to perform the finite element analysis of the waveguide structure by placing 
the structure inside a computational window, which is large enough to guarantee the field vanishing at the boundary. 
Here, we consider an optical wavelength of 1 fim and the optical permittivities of the waveguide are eq — —50 — 3.85i 
and £a = 2, corresponding to gold and polymer^^. The dispersion and the field orientation of the possible modes for 
cylindrical and square waveguides are presented in Fig. [5] and Fig. [3] respectively. As shown in the inset of Fig. [51 
these modes can be presented by two indices, where the first index denotes the number of the angular moment, m, 
and the second index describes the polarization degenerate mode with the same m. For example, if the E'^'^ denotes 
the mode with angular moment of m, then E"^'^ denotes the corresponding degenerated mode, the field distribution 
of which is rotated by 9 along the z axis compared with where 9 = 27r/2m. As pointed out by Takahara et 

ali^, the fundamental mode E^^'^^^ does not have a cutoff size of the radius, which is confirmed from the dispersion 
relation in Fig. [21 The modes supported by the metallic nanowire preserve the cylindrical symmetry of the waveguide. 
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Due to the constraints from the boundary conditions, only TM modes exist. For the square surface plasmon-polariton 
waveguides, the fundamental modes, which were studied by Jung et ali^, can be labeled in terms of two indices, which 
denote the number of sign changes in the dominant component of the electric field along the x and y axes respectively. 
Both plasmonic waveguides support one fundamental mode {E^'^) without any cutoff size of the metal core, and the 
corresponding propagation constants increase when the size of the metal core is further shrunk, which slows down the 
propagating plasmonic mode. Such geometric slowing down enhances the local density of the states and the coupling 
efficiency to the nearby quantum emitter. In the following calculations, the size of metal core is restricted below the 
cutoff size of higher order modes so that only a single mode is supported. 

The electric-field dyadic Green's function for a specific guided plasmonic mode is constructed from the numerical 
calculation of the electric field. In the following part we will explain how to construct the electric-field dyadic Green's 
function for one guided plasmonic mode^S. 

The electric dyadic Green function G(f , uj) is defined by 

[V X V X ~kle{f)]G{f, u) = 16{f - f'), 

where / is the unit dyad. Rigorously speaking, the operator defined by L = [V x V x — fcoe(f)] does not have a set 
of complete and orthogonal eigenmodes due to its non Hermitian character. Without loss of generality, we adopt 
biorthogonality in the present paper to form a complete set of "orthogonal" modes of the waveguides initially, and 
then we will end up with an approximation from the power orthogonality for our plasmonic waveguides. Suppose 
that En are a set of eigensolutions defined by L, the biorthogonal modes E^ are defined as the eigensolutions of the 
adjoint operator denoted by U , which is obtained from the operator L by replacing e{f) with its complex conjugate. 
The biorthogonality condition is given by 

e(r)En{r)[El{f)]*d\ = S^mN^, (1) 
with the completeness relation ^ ^('')-^"(^[-^..('' )1 — JS{f — f'). From the biorthogonal completeness relation, the 



dyadic Green function G{r^uj) can be constructed from the eigenfunction expansion as follow! 
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G{f,f') = GGT{f,f') + GGL{f,f') 

_ ^ E^(f)-[El(f')]' ^ V0^(f)-[v4(f')]-' (2) 

— Ar„A„ MiTfc^ 

n n 

where the generalized transverse part of the dyadic Green's function, Ggt, is constructed from the complete set of 
transverse eigenfunction i?„ (f) given by, 

-V X V X En{f) + kle{f)E^{f) = Ke{f)En{f), 
V • Hr)En{f)] = 0, 

with the eigenvalue A„. The longitudinal or quasistatic part Ggl is constructed from longitudinal eigenfunction 
that can be found from a complete set of scalar eigenmodes </)„ (f) satisfying 

V- [e(f)V0„(f)] =a„<^„(f) (4) 

with the biorthogonality relation, J e(f)V0„(f) • [V0jj(f)]*(i'^r = SnmMn- Since we are studying the guided plasmonic 
mode, which describes the field solution in the absence of electric charge (V • [e(r)i?„(f)] = 0), the longitudinal 
component will vanish in the following calculations. 

By applying the principle of constructing the electric-field dyadic Green's function to the case of a plasmonic 
waveguide, we find the contribution to the dyadic Green's function from the plasmonic modes as 

^ — oo 

where a = {p,—j3}, and the normalization factor N is given by S{P — (3')Sppi N = J e{f)E{f)[E'' {f)]*(Pr = 
2'k5{(3 — P')Spp' J e{f)Ea{r)[El, {f)]*dxdy, which can be further simplified as = 27r / e{f)Ea{r)[El, {f)]*dxdy, where 
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a denotes {p , — /3 }. For one plasmonic mode, the expression (5) is evaluated in closed form by the method of contour 
integration as the integrand decays to zero at infinity in the upper and the lower /3 plane, 



UjNVg ' 

where Vg is the group velocity, defined by Vg = duj/dp. The corresponding density of states 
for one plasmonic mode can be calculated from the dyadic Green function according to Novotny^, 
Pp(ro, ujq) = 6uj[n^ ■ Im{G(ro, rg, i^o)} ■ ftf^]/ {ttc^), where is the unit vector of the dipole moment. If the dipole emit- 
ter is oriented along the X axis, the density of states for the plasmonic mode is given by Ppi{f, ui) = 6\Ea^x{r)\'^ / (Nvg). 
The spontaneous emission decay rate into the plasmonic mode can be calculated by -fpi = j^\fi\'^ppi{f,uj). Normal- 
ized by the spontaneous emission decay rate in the vacuum, the emission enhancement due to the plasmonic excitation 
is 



- (7) 
70 (^iNvg 

Eq. ([7]) gives a general expression of the spontaneous emission decay rate into a guided mode, supported by a lossy 
or lossless waveguide. In dielectric waveguides, losses are generally small, and the biorthogonal modes E^^ can 
approximately be replaced by the orthogonal mode Em- Such an approximation is also valid for our plasmonic 
waveguide, where the imaginary part of the propagation constant for the fundamental mode is around 1% of the real 
part. According to Snyder—, the group velocity can be calculated by Vg = {Ex H*) ■ zdAj eoe{r)\E{r)\'^dA, 
where Ago denotes integration over the transverse plane. By applying the power orthogonal approximation and the 
explicit form of the group velocity to Eq. ([7]), we can reach the following expression for the plasmonic decay rate of 
the fundamental mode, 

7p[ ^ 37rc£o^ao,x(r)^;„^j^(f) 

70 ~ k^lA {ExH*)-zdA- ^ ' 



B. Total decay rate 

As described in the previous subsection, the well defined field components in the transverse plane of the waveguide 
give the possibility of constructing the dyadic Green's function numerically. The reason is that the field is concentrated 
around the metallic core and is decaying to zero on the borders when the modeling domain is reasonably large. Hence, 
the perfect electric conductor boundary condition is implemented to truncate the 2D modeling domain. However, 
for the radiation modes, the field components in the transverse plane of the waveguide do not vanish no matter how 
large the modeling domain is. Hence, it is extremely difficult to construct the dyadic Green's function numerically 
for the radiation modes in a similar way as for the guided mode. Therefore, we implement a 3D model to include 
the radiation modes, as well as the nonradiative contributions, by solving the wave equation with a harmonic (time 
dependent) source term, 

[V X —V X -kle{f)]E{f, Lu) + jLunoJiui) = 0. (9) 

fir 

If we introduce a test function F(f, w), we can construct the functional corresponding to the wave equation in the 
following way^, 

^ = JJJ ["^ X —V X ~k'^e{f)]E{f,uj) ■ F*{f,Lu)dV + j j j jujfioJ{Lo) ■ F*{r,uj)dV 

V V 

—V X E{f,uj) - V X F*{f,uj)dV ^ j j j k^e{f)E{f,Lu)-F*{f,u)dV + JJJ jujfxj{uj) ■ F*{f,uj)dV (^g) 

V V V 

+ (^F* {f,uj)-[—nxV X E{f, uj)]ds , 
av 
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FIG. 4: A single quantum emitter coupled to a metallic nanowire. The grey transparent region represents the perfectly matched 
layers, the mode matching boundary condition is applied on the top and the bottom of the structure. The quantum emitter is 
implemented by an electric line current. 

where dV denotes the surface that encloses the volume V , and n denotes the outward unit normal vector to the surface 
of the modeling domain,. This is the variational formulation of the wave equation, which is required to hold for all 
the test functions. Eq. (jlOp enables us to formulate the finite element solution for such a boundary-value problem 
by employing the standard finite element solution procedures, including discretization and factorization of a sparse 
matrix"^^. The boundary- value problem defined by Eq. (jlOp was solved by utilizing a commercial software package, 
COMSOL Multiphysics^^. 

It is crucial to truncate the computational domain properly. As shown in Fig. 01 we have two strategies to truncate 
the modeling domain: I) In the X-Y plane, the computation domain is truncated by the perfectly matched layers 
with thickness of half a wavelength in vacuum. II) Along the Z-axis, the computation domain is terminated by mode 
matching boundary conditions, which will induce a certain amount of reflection from the radiation mode and the 
higher order plasmonic modes if they exist. Essentially, the mode matching boundary is an absorbing wall, which 
behaves as a sink of electromagnetic waves. There are different options of realizing the mode matching boundary to 
absorb a single mode, depending on whether the absorbed mode is TE, TM or a hybrid mode. For a pure TM or TE 
mode, it can be matched by simply applying the conditions, 



—n X V X E{f,uj) = ° y \ 

— n X V X E{f, uj) — jPn x — n x Et{f, uj), 



TM; 
TE- 



(11a) 
(lib) 



on the boundary, where /3 is the propagation constant, and Et{f,uj) is the tangential components of the dependent 
variable E{f,u}) on the boundaries in the numerical model. The mode matching boundary condition for the hybrid 
mode can be implemented as 



1 

Mr 



n X V X E(r,uj) = ^juj^Qn x iJo, 



(12) 



where E{f^ lo) is the dependent variable solved in the 3D model, and Hq denotes the matched mode that is applied. In 
our model, Hq corresponds to the fundamental hybrid mode supported by the plasmonic waveguide. It is calculated 
from the 2D eigenvalue problem, and is given by 



^2d 



{Hox, Hoy, Hoz). 



(13) 



Here, P2d, and f3 are the time averaged power flow, the magnetic field, and the propagation constant, respectively 
calculated from the 2D model, while Pq denotes the normalization factor of the power emission in the 3D model, 
and Lq represents the half length of the 3D model. Due to the losses of the metals, the magnitude of the magnetic 
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TABLE I: The relation of the 6 field components for the fundamental hybrid mode 



Description 


Relation 


Tangential electric field, 


ElD,t = £^in X H2D,t)s + ^{\7t X 


s G [x,y] 




Normal electric field 


E5d = ■ {Vt X H2D,t) 


Tangential magnetic field, 


^2D,t — J^2d,m 


s G [x,y] 




Normal magnetic field 


H2D = jH2d,m 
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FIG. 5: Length dependence study of the total decay rate for the metallic nanowire. The radius of the metallic nanowire is 20 
nm, the distance of emitter to the wire edge is 30nm. 



field is a complex number. In order to guarantee that the phase of at the position of the emitter is zero when 
the emitter is oriented horizontally, the extra phase = arctan( '™°f(^J'x ° ) needs to be compensated, i.e., Hq = 

\/ '^P2d° H2d&~''^^'"°^'^'' ■ In the 2D eigenvalue calculations, there are 6 components involved for the hybrid fundamental 

model, the relations of which are tabulated in Table HI The magnitudes of the magnetic field, {^H^j^ Hl^j^ ^, ^ , 
are the dependent variables, which are calculated directly from the 2D numerical model. 

The total decay rate, jtotai, is extracted from the total power dissipation of the current source coupled to the 
nearby metallic waveguide, jtotai/jo = Ptotai/Po, where Ptotai = 1/2 Jy Re{J* ■ Etotai)dV is the power dissipation of 
the current source coupled to the metallic waveguide, and Pq — 1/2 jy Re{J* ■ Eo)dV is the emitted power by the 
same current source in vacuum. Pq is a normalization factor, which is also used to normalize the power flow on the 
boundaries in Eq. (|13p . As demonstrated in Fig. [H the field is generated by the current source, namely, the dipole 
emitter, which is implemented by a small electric line current. In our model, the dipole is oriented horizontally. For 
an electric current source with finite size of Z <C Aq), and linear distribution of current /q, the dipole moment of the 
sourceS is, /i = jIqI/uj. In order to avoid higher order multipole moments, the size of the current source should be 
restricted below a certain value. Our numerical test shows that the variation of the total power dissipation from the 
size dependence of the emitter is negligible when the size of emitter is shorter than 2 nm. 

In order to check the validity of the mode matching boundary condition we studied the length dependence of the 
total decay rate for two different plasmonic waveguides. The length dependence of the total decay rate ptotai for the 
metallic nanowire is shown in Fig. [S] The fundamental mode supported by the metallic nanowire is TM, hence the 
mode matching boundary condition defined by Eq. (jllap is implemented. As can be seen from Fig. [SJ the variations 
in the total decay rate are reduced by increasing Lq, and the damped oscillation of the total decay rate with Lq 
indicates a certain amount of reflection from radiation modes, which is confirmed by the period of the oscillation ( 
equal to the wavelength in the media with e — 2). We also see that the variation of the total decay rate due to the 
length dependence is below ±2.5% due to the dominating excitation of the plasmonic mode when Lq is larger than 1 
/iTO. Basically, the accuracy of ptotai Mo relies on the length of the plasmonic waveguide, accordingly we estimate the 
relative error on the computed data is ±2.5% in the following calculations for the metallic nanowire. 

Regarding the square plasmonic waveguide, the condition defined by Eq. (|13p is applied on the boundary to absorb 
the hybrid mode supported by the waveguide, where Hq is the magnetic field for the matched field. As shown in Fig. [51 
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FIG. 6: (a) Length dependence of the total decay rate for the square plasmonic waveguide. The side length is 30nm, the 
distance of the emitter to the edge of the square metal core is 20 nm. (a) Length dependence study with damped oscillations, 
(b) Length dependence for the points from (a) (marked ellipses) where where real{e~-'^^'^°'^'^^) = holds approximately, (c) 
Illustration of the reflection of the normal magnetic field of the fundamental hybrid in the 3D model, r and 9 are the reflection 
coefflcient and phase shift respectively. 



there is also a damped oscillation of the total decay rate with the length of the computation domain, and the tendency 
of achieving higher accuracy for 740*0; when Lq is lengthened, which is similar to the length dependence study of the 
total decay rate for the nanowire. Nevertheless, there are two distinctions between the two plots: I) The variation 
of the total decay rate for the square plasmonic waveguide is much larger than that for the metallic nanowire; II) 
The variation of the total decay rate for the square plasmonic waveguide with Lq primarily stems from the reflection 
of two different modes, which are indicated by two different periods in the damped oscillation. The reflection of the 
fundamental mode, which is supposed to be absorbed at the boundaries, is responsible for the oscillation with the 
period of 200nm, the other oscillation with the period of 370 nm results from the reflection of a quasi guided mode, 
denoted by Egg . The explanation is the following, the boundary condition defined by Eq. (jllap can completely absorb 
the matched pure TM mode, while it is not true for the boundary condition defined by Eq. (|13p for the hybrid mode, 
and a significant reflection from a quasi guided mode also exists for the square plasmonic waveguide. For the hybrid 
mode the last term in Eq. (|10p relies not only on the tangential components of the electric (magnetic) field, but also 
on the normal component of the electric (magnetic) field, which is intrinsically lost on the boundary in the vector 
element formulation of the 3D numerical modeP'^. Our interpretation is that, even though the normal component 
of the electric field can be included on the boundaries by Eq. (jl3p . the normal component of the magnetic field is 
essentially missing in the 3D numerical model with the square plasmonic waveguide, resulting in the reflections in 
our vector element formulation'^"''"^^. However, in Fig.[Sfa), it appears that the points for which real{e~^^^^°~^'^'^) — 
holds approximately converge quickly with minimum impact of the reflection from the fundamental hybrid mode. The 
mode Egg, with effective wavelength of 740.07 nm, is characterized by the material properties of the waveguide and is 
rather insensitive to the size of the metallic core. Compared with other quasi guided modes or radiation modes, the 
mode Egg has a relatively significant contribution to the jtotah the normalized spontaneous emission rate is 0.107. 
Since no extra effort is made to prevent the reflections from any components of the mode Egg, it is understandable 
that the induced reflections give rise to several peaks in Fig. El^a). 
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The normal component of the magnetic field of the fundamental mode in the 3D model can be obtained by a 2D 
eigenvalue calculation, H^^i = \J '^p^° H2^e~^^^''^'^^ , where / is the distance from the observation plane to the emitter. 
Similarly, the reflected normal component of the magnetic field at the position of the emitter can be obtained by 
taking into account the phase shift due to propagation and reflection, H^q — t\J '^'p^" ffn^Q-j{^0Lo+4>+s) ^ shown 
in Fig. [6] (c). The reflected normal component of the magnetic field will "generate" a perturbation term EJ^ to the 
original component, the real part of which is integrated to calculate the total power dissipation. According to 
Table HI the reflected term i?^ from the fundamental hybrid mode is given by 

= "i^^* ^ -(-Y^i?k™e-^-(^'^^°+^+^))).. (14) 

The real part of E^ can be zero when Lq is appropriately chosen, therefore, the obtained total decay rates are expected 
to approach the true value more closely due to the vanishing contribution of E^ to the total decay rate. In Fig. [61(a), 
at the points with marked ellipses, the half model length Lq fits the requirement {real{E!^) = 0), and we also found 
that the phase shift 9 is required approximately to be 7r/2. Further examining the phase shift 9 involves technical 
details regarding the implementation of the vector element formulation of the finite element method, which is beyond 
the scope of the present paper, and we refer to the references22i2ii^i2&. From Fig. [HKb), we estimate the relative error 
on the computed data for the square plasmonic waveguide to be ±2%, when Lq is larger than 1 /xm. 

III. RESULT AND DISCUSSION 



In this section, the numerical method is applied to the two cases, the metallic nanowire and the square plasmonic 
waveguide. For the metallic nanowire, we compared our numerical calculations with the quasistatic approximation, 
which was studied by Chang et al2i. As can be seen in Fig. [7l our numerical results agree well with the quastatic 
approximation when the radius of the metallic nanowire is less than 20 nm, while it is 5-10 times larger when the 
radius is 100 nm. The deviation between the two methods increases when the radius becomes larger. The results of 
the comparison can be understood if one realizes that the quasistatic approximation is valid only when the radius is 
much smaller than the wavelength. For the large wires, which has also been pointed out by Chang et ali^, the full 
electrodynamic solutions predict significantly larger values of "fpi/"fo and the /? factor. The quasistatic approximation 
assumes that the magnetic field vanishes, thus the obtained solution of the electric field simply behaves as a static 
field. It is the same for the plasmonic mode, the penetration length of which in the dielectric medium is considerable 
shorter than that obtained from the full electrodynamic solutions, and therefore the quasistatic calculation predicts 
a lower coupling efficiency. In summary, our numerical calculation is consistent with the quasistatic approximation 
for the nanowire radii approaching 0, and the values from finite element simulation are generally larger than those 
obtained from the quasistatic approximation when the radius is beyond 20 nm. 

We also studied the coupling of the quantum emitter with the square plasmonic waveguide. As shown in the inset 
in Fig. [51 the quantum emitter is oriented along the X axis, and the distance dependence of the plasmonic decay rates 
and spontaneous emission (3 factors is calculated as function of emitter position along the X axis. With optimized 
side length of the waveguide and distance of the emitter to the edge of the waveguide, the /3 factor can reach 80%. 

IV. CONCLUSION 

In conclusion, we developed a self-consistent model to study the spontaneous emission of a quantum emitter at 
nanoscale proximity to a plasmonic waveguide using the finite element method. The dyadic Green function of the 
guided modes supported by the plasmonic waveguide can be constructed numerically from the eigenmode analysis, and 
subsequently the normalized decay rate into the plasmonic channel can be extracted. The 3D finite element model is 
also implemented to calculate the total decay rate, including the radiative decay rate, nonradiative decay rate, and the 
plasmonic decay rate. In the 3D model, it is assumed that only one guided plasmonic mode is dominatingly excited, 
which is normally true when the size of the cross section of the plasmonic waveguide is below 100 nm. Under such 
condition, the spontaneous emission p factor is calculated. We compared our numerical approach with the quasistatic 
approximation for the gold nanowire. We observe agreement with the quasistatic approximation for radii below 20 nm, 
where the quasistatic approximation is valid. For larger radii the FE simulation predicts approximately 5 times larger 
values. This is reasonable since the numerical model takes into account wave propagation, whereas the quasistatic 
approximation calculates the static field. We also applied our numerical model to calculate the spontaneous emission 
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FIG. 7: Comparison of FEM simulated results based on the dyadic Green function with the quasistatic approximation for the 
metallic nanowire. 
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FIG. 8: Distance dependence of the plasmonic decay rates and spontaneous emission fi factors for the square plasmonic 
waveguide. 



of a quantum emitter coupled to a square plasmonic waveguide. The numerical calculations shows that spontaneous 
emission /3 factor up to 80% can be achieved for a horizontal dipole emitter, when the distance and the side length 
are optimized. 
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